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Abstract 


For scattering calculations involving nonspherical particles such as ice crystals, 
we show that the transverse wave condition is not applicable to the refracted 
electromagnetic wave in the context of geometric optics when absorption is involved. 
Either the TM wave condition (i.e., where the magnetic field of the refracted wave is 
transverse with respect to the wave direction) or the TE wave condition (i.e., where the 
electric field is transverse with respect to the propagating direction of the wave) may be 
assumed for the refracted wave in an absorbing medium to locally satisfy the 
electromagnetic boundary condition in the ray tracing calculation. The wave mode 
assumed for the refracted wave affects both the reflection and refraction coefficients. As 
a result, a nonunique solution for these coefficients is derived from the electromagnetic 
boundary condition. In this study we have identified the appropriate solution for the 
Fresnel reflection/refraction coefficients in light scattering calculation based on the ray 
tracing technique. We present the 3x2 refraction or transmission matrix that completely 
accounts for the inhomogeneity of the refracted wave in an absorbing medium. Using the 
Fresnel coefficients for an absorbing medium, we derive an asymptotic solution in an 
analytical format for the scattering properties of a general polyhedral particle. Numerical 
results are presented for hexagonal plates and columns with both preferred and random 
orientations. The asymptotic theory can produce reasonable accuracy in the phase 
function calculations in the infrared window region (wavelengths near 10 pm) if the 
particle size (in diameter) is on the order of 40 pm or larger. However, since strong 
absorption is assumed in the computation of the single-scattering albedo in the 
asymptotic theory, the single scattering albedo does not change with variation of the 
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particle size. As a result, the asymptotic theory can lead to substantial errors in the 
computation of single-scattering albedo for small and moderate particle sizes. However, 
from comparison of the asymptotic results with the FDTD solution, it is expected that a 
convergence between the FDTD results and the asymptotic theory results can be reached 
when the particle size approaches 200 pm. We show that the phase function at side- 
scattering and backscattering angles is insensitive to particle shape if the random 
orientation condition is assumed. However, if preferred orientations are assumed for 
particles, the phase function has a strong dependence on scattering azimuthal angle. The 
single-scattering albedo also shows very strong dependence on the inclination angle of 
incident radiation with respect to the rotating axis for the preferred particle orientations. 
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1 . 


Introduction 


The optical properties of dielectric particles such as aerosols and ice crystals in 
the atmosphere are fundamental to a number of disciplines including atmospheric 
radiation transfer and air-borne or satellite-borne remote sensing applications. In a recent 
book 1 , Mishchenko et al. have comprehensively reviewed various methods that have been 
developed to solve the scattering and absorption properties of particles for a variety of 
geometric morphologies and electromagnetic characteristics of the scatterers. Accurate 
numerical methods such as the discrete dipole approximation (DDA) method 2-3 and the 
finite-difference time domain (FDTD) technique 4,5 are applicable in practice only to size 
parameters smaller than 20 because the computational requirements increase quickly with 
size parameter. Although analytical solutions are available for some particle shapes such 
as spheres 6 and spheroids 7-8 , the corresponding numerical computations are usually very 
challenging. For example, the computation of the optical properties of spheroids based on 
an analytical solution in terms of a series of special functions may not be stable when the 
size parameter is larger than approximately 30-40. It is worth noting that the recent 
developments on the T-matrix method 9-10 allow the exact solution for the optical 
properties of spheroids and finite circular cylinders for size parameters up to 200 11 . There 
is no single method to solve for the optical properties of nonspherical particles across the 
entire size parameter spectrum. 

It is common to use the ray-tracing technique 12 ' 15 for nonspherical particles that 
have sizes much larger than the incident wavelength. Previous applications of the ray- 
tracing technique have been focused mainly on scattering computations for the visible 
and near-infrared wavelengths for which the imaginary part of the refractive index is 
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small. For these weakly absorptive spectra, the conventional Snell’s law and Fresnel 
formulae can be used to account for the ray reflection and refraction at particle surface. 
However, the optical properties of large dielectric particles at wavelengths with strong 
absorption are fundamental to many practical applications. For example, the scattering 
and absorption properties of ice crystals at infrared wavelengths are required in the 
retrieval of cirrus clouds using satellite radiance data such as that from the Advanced 
Very High Resolution Radiometer (AVHRR) 16 and Moderate Resolution Imaging 
Spectroradiometer (MODIS) 17 imagers. Under cirrus cloudy condition, interpretation of 
the High resolution Interferometer Sounder (HIS) 18 hyperspectrum in the atmospheric 
window region (8-12 pm) also requires the optical properties of ice crystals. 

When ray optics is applied to the infrared spectrum in regions of strong 
absorption, the refracted wave inside the particle is inhomogeneous in the sense that the 
planes of constant phase and constant amplitude are not parallel 19 . Current applications of 
the ray tracing technique assume that planes of constant phase and constant amplitude are 
parallel. The inhomogeneity associated with the reflected wave substantially complicates 
the behaviors of rays at particle surface and inside the medium. Stratton 19 and Bom and 
Wolf 20 have discussed this inhomogeneity and its effect on the Fresnel reflection 
coefficients associated with complex refractive index. In particular, Stratton 19 (section 
9.9) has solved the problem of the reflection of a plane electromagnetic wave by a 
conducting surface. Recently, Zhang and Xu 21 and Yang and Liou 22 addressed the effect 
of the inhomogeneous wave on ray-tracing calculation. Zhang and Xu 21 have expressed 
the Fresnel coefficients without an explicit use of complex refractive index. However, 
they did not account for the effect of the inhomogeneity on wave attenuation that is the 
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most important physical process involved in the ray propagation within an absorbing 
medium. Yang and Liou 22 have discussed the impact of the inhomogeneous effect on ray 
propagation by introducing an effective refractive index, but they employed the high- 
frequency approximation 20 in the calculation of reflection and refraction coefficients. 
Since the effect of inhomogeneity on ray tracing calculation has not been solved 
completely, it is necessary to further study this problem. The intent of this study is to 
develop an analytical asymptotic method that comprehensively accounts for the 
inhomogeneity effect when deriving the optical properties (such as the single-scattering 
albedo and phase function) of strongly absorbing particles with sizes much greater than 
the incident wavelength. 

To compute the optical properties of dielectric particles with strong absorption we 
first derive the reflection and refraction coefficients for the electric field in terms of an 
apparent refractive index without using the high-frequency approximation. The reflection 
and refraction coefficients incorporate the inhomogeneous effect on the direction of ray 
propagation and the Fresnel coefficients. Subsequently, we derive the asymptotic solution 
expressed in an analytical form for large particles with strong absorption. In Section 2 we 
discuss the various wave modes of an inhomogeneous wave refracted into a particle when 
absorption is involved. We show that the electric or magnetic field may not be 
perpendicular to the direction of ray propagation because of the effect of strong 
absorption. The Fresnel coefficients are reformulated, in particular, the transmission 
coefficient for the component of the refracted field that is parallel to ray propagation. In 
Section 3, we present the analytical solution for the optical properties of large polyhedral 
particles with strong absorption. The solution is derived on the basis of the 
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electromagnetic relationship between the near field and far field. In Section 3, we also 
present numerical results of the asymptotic solution for hexagonal ice crystals. Finally, 
conclusions of this study are given in Section 4. 


2. Physical Basis for Ray Tracing in Absorptive Case 

In this section, it is shown that the transverse wave condition cannot be applied to 
both the electric and magnetic fields associated with a ray refracted into an absorbing 
medium. In addition, we show the wave mode subject to an electric boundary condition is 
not unique in the context of geometric optics when absorption is involved. The behavior 
of the refracted wave not only affects the Fresnel refraction coefficients but also the 
Fresnel reflection coefficients. In the following discussion, it is shown that there are four 
combinations possible for the reflection/refraction coefficients that depend on wave 
modes. While each combination of the reflection and refraction coefficients is a 
mathematically rational solution, only one set of coefficients should be used in the ray 
tracing solution. The choice of the proper wave mode for the refracted wave in ray 
tracing calculation is made through the comparison of the geometric optics solution to the 
analytical Mie result for the scattering properties of a very large (size parameter on the 
order of 1000) sphere with strong absorption. 

Let us consider the application of the geometric optics method to solve for the 
scattering properties of a particle that is large and absorptive. In this case, the rays that 
are refracted into the particle are largely absorbed. Thus, we need only consider the first 
order reflected/refracted rays because of the small amount of energy carried by the rays 
that undergo internal reflections. To examine the first order reflection/refraction 
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characteristics and the corresponding polarization features of the electric field, we begin 
with the phase variation for the field associated with a complex refractive index. Let three 
unit vectors, e n e r , and e t , indicate the propagating directions for the incident, reflected, 
and refracted rays. To systematically present this study, it is necessary to recapitulate 
some results obtained by Yang and Liou 22 . Following these authors, the fields associated 
with these rays can be written as follows: 


E i (r) = E io exp(ike i r). 

(la) 

E r (f) = E ro exp(ike r r), 

(lb) 

Et (0 = E t0 exp[ik(N r e, + iN n n) ■ r]. 

(lc) 


where E i o , E r o , and E t0 are the amplitude vectors for the fields; k = 2n/X in which X 

is the wavelength in vacuum; n is a unit vector that is normal to particle surface and 
points into the particle. The two parameters N r and N n in Eq.(lc) are given by 

N r = l~ ll2 {m r 2 - m 2 +sin 2 ^,. + [(m 2 - m 2 -sin 2 £,.) 2 + 4m 2 m 2 ]' /2 }, and (Id) 

N n = 2" l/2 {-(m 2 - m 2 - sin 2 Q + [{m 2 - m 2 - sin 2 Q 2 + 4m 2 m 2 ] 1/2 } , (le) 

in which m r and m j are the real and imaginary part of the refractive index, and is the 
incident angle given by ^ = sin -1 [l-(e ( . f -n) 2 ] 1 ' 2 . Note that N„ (denoted as N* in Ref. 22) 
is not presented in the previous study. Physically, N r is the parameter that determines the 
phase and N n is the parameter that determines the absorption. This will be illustrated in 
the following discussion. 

From Eqs.(la)-(lc) and the continuity of wave phase at the medium interface, it 
follows that Snell’s law for an absorbing medium can be mathematically expressed in the 
form of 
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A A . A/A AvA 

= e,. + 2(e i ■ n)n, 

(2a) 

e, = e,. / N r + (cos£, - cos£,. / N r )h, 

(2b) 


where the refractive angle is given by = sin '[sin(i^)/ iVJ. Equation (2a) means that 
the incident and reflected rays are spatially symmetric about the direction normal to 
particle surface at the incident point, which is the standard Snell’s law for reflection 
direction. However, Eq. (2b) indicates that the refracted direction is not determined by 
m r but by yV r . For this reason, N r is referred to as the real part of apparent refractive 
index, which depends on the incident configuration and the dielectric constant of the 
medium for the determination of the direction and phase for refracted wave. The 
attenuation of the refracted field due to the absorption depends on the path of 
observation. If the path is along n, the attenuation is determined by N n , as is evident 

from Eq.(lc). However, if we trace the refracted wave along the direction of the refracted 
ray, i.e., along position vector r - le, in which / is the penetration depth of the ray, the 

electric field is given by 

E t (le t ) = E t0 exp [~kN n (h ■ e t )/] exp (ikN r l) . (3) 

Obviously, the factor NJh-e t ) is the imaginary part of the apparent refractive index that 

governs the attenuation of the electric field associated with the rays. We denote this 
factor as N t . Thus, for a ray impinging on an absorbing medium, the apparent refractive 
index in the complex format is (N r + /A',). It can be proven that 

Nj = N„(n ■ e t ) = N n cos£, = m r m i / N r . (4) 

Figure 1 shows the real and imaginary parts of the inherent and apparent 
refractive index of ice at infrared wavelengths 1 1 and 12 pm for which substantial 
absorption exists. The optical properties of ice at these two wavelengths are important 
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for many applications in satellite remote sensing 17 . Note that for ice particles, 
wavelengths near 1 1 pm are in the Christiansen band where absorption dominates the 
extinction 23 . The inherent refractive index is determined uniquely by the dielectric 
characteristics of the medium, i.e., permittivity, which does not depend on the incident 
configuration of a ray. As is shown in Fig. 1 , the real part of the apparent refractive index 
increases with the increase of incident angle whereas the imaginary part of the refractive 
index decreases with the increase of the incident angle. Obviously, in terms of wave 
attenuation per unit length along propagating direction of a ray, the particle is less 
absorptive at a large incident angle. For incident angles larger than 10°, the 
inhomogeneous effect due to strong absorption cannot be ignored. 

It can be shown that various wave modes may satisfy the electromagnetic 
boundary condition associated with the reflection and refraction at the interface of two 

media, one of which is strongly absorptive. First, let P be a unit vector perpendicular to 
the incident plane. The unit vectors that are parallel to the incident plane and 
perpendicular to ray propagation directions for the incident, reflected, and refracted 
waves are given by 

a,. =e,. xj3, a r =e r xp, and a, =e,xp. (5) 

It should be pointed out that pd f r ,e f t constitute right-handed coordinate systems. The 

incident and reflected waves are within non-absorbing medium. As a result of the 
homogeneous wave properties of the incident and reflected rays, transverse wave 
condition can be imposed on both electric and magnetic fields for these rays (hereafter, 
referred to as TEM mode). Thus, the corresponding amplitude vectors of the electric 
fields can be decomposed in the form of 
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(6a) 


4 = 4>,+ 4,#P. 

E ro = E ro,o^r + ^ro, fiP * 


(6b) 


The Maxwell equations in time-independent form are given by 


E (r ) = V X H(r)/(-iks), 

(7a) 

ff(f) = V x E(r) / ik , 

(7b) 


where £ is permittivity given by £ = (m r +im i ) 2 . We have chosen the time harmonic 
factor as exp (ikct) in which c is the speed of light in vacuum. If this factor is selected as 
Q\p(-ikct), the sign of the imaginary part of refractive index is negative 1,24 . Using 
Eq.(7b), we can obtain the amplitude vectors of magnetic fields associated with the 
incident and reflected rays as follows: 


H to = e, x (E, o a a, + E w $) = E i0 fi a, - E io a P, 


K =e r * ( E ro,A + E ro# P) = E roA ~ E ro,a P‘ 


(8a) 

(8b) 


For the refractive wave in an absorbing medium, a transverse wave condition on both the 
electric and magnetic fields cannot be imposed because of the inhomogeneity effect, that 
is, the TEM wave mode is not valid in this case. 

In the following discussion, we first examine the reflecting characteristics of an 
absorbing medium for various wave modes. To derive the reflection coefficient, the 
electric field may be assumed to be transverse with respect to the propagating direction of 
the ray (hereafter, this mode is referred to as TE mode) so that the refracted electric field 
can be expressed in the form of 


E ,o,,A, + E <o, pi 


p. 


( 9 ) 
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From Eqs. (lc), (4), (7b), and (9) and the relationship of n = cosine, -sin£ ( a,, the 
refracted magnetic field in TE mode can be obtained as follows: 


HJf) = (N r e, + iN.fi) x (£„>, + EJft 


= W + - W + + W. sin?, • 


( 10 ) 


From Eq.(10), it is evident that the refracted magnetic field has a nonzero component 
along the unit vector e, (ray direction) when N n is not zero, i.e., when absorption is 

involved. At the interface of two media, the electromagnetic boundary condition requires 
that the tangential components of electric and magnetic fields are continuous. Thus, we 
have the following relationships: 

E iofi + E ro £ = E t0 #’ ( 1 1 a ) 

E iofi cos;,. - E ro? cos;,. = [(AT + /AQcos;, + iN n sin 2 (1 lb) 

= W+w,)£„, O'<0 

COS?, - cos?, = cos ?,£„,,. (lid) 

The preceding equations may be solved to obtain the reflection coefficients for the two 
polarization configurations: 

R „ /F _ cos;,. - [(N r + /A^Qcos;, + iN „ sin 2 ;,] 

TE,L ro ,p (o,p cos ; . + [(N r + iN t ) cos + iN n sin 2 ;, ] 


_ cos;,.-(iV r cos;, + /iV n ) 
cos^.+^cos;, +iN„y 


n _ F IP _ (A r +/A,.)cos;,.-cos;, 
TEJ! ro.O. ro# (AT+/JV.) cos;, + cos;, 


(12a) 

(12b) 
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Similarly, if we assume that the magnetic field of the refracted wave is transverse with 
respect to the ray direction (hereafter, this mode is referred to as TM mode), the 
reflection coefficients are then given by: 


= (N L ±iN l ) cosC-ecosC ^ 

m± (N r + iNJcosC,! + ecos£ ( 


R _ ecos^,-(N r cos^,+iN m ) 
™’" £ cos - (N r cos + iN n ) 


(13a) 


(13b) 


Therefore, considering the two wave modes and two polarization configurations, four 
combinations are possible for the reflection coefficients: (R m// ,R m ^), {Rte,h^te,l)^ 
(i? r£ // ,i? rA/1 ), and (R mj ,,R TE ' ± )- Each of these four combinations is a mathematically 
rational solution derived from the electromagnetic boundary condition when an 
absorptive medium is involved. However, only one among the four solutions should be 
computationally correct in the ray tracing calculation regarding the process of scattering 
of electromagnetic wave by a dielectric particle with absorption. It should be pointed out 
that the Maxwell equations subject to appropriate boundary conditions always have a 
unique solution, as is in the case of Mie theory. The explanation for the multiple solutions 
for the reflection coefficients is that electromagnetic boundary condition is locally 
imposed in the context of geometric optics, which allows the treatment of the 
inhomogeneous refracted wave as a refracted ray, thereby providing the basis for a ray- 
tracing procedure. Accordingly, we are looking for the pair of reflection coefficients that 
provide the best numerical agreement with the exact solution in scattering computation. 

To identify the computationally correct pair of reflection coefficients among the 
preceding four combinations, we consider the scattering geometry in Fig. 2 for a strongly 
absorptive sphere with a size much larger than incident wavelength. In this case, the 
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scattered field is composed primarily of the contributions from diffraction and external 
reflection. The incident energy associated with the area element ds on the projected area 
is scattered in angular element sin QdQdq in which 0 is scattering angle. The scattered 
intensity at a distance r that is far away from the particle is given by 

\R\ 2 a 2 


/ = / 


l*l 2 

ds 

r 2 

sin QdQdty 


= / 


4r 


,2 ’ 


(14) 


where a is the radius of the sphere; I 0 and I s are incident and scattered intensity, 
respectively; R is reflection coefficient without a specification of wave mode and 
polarization configuration. To derive Eq.(14), we have used two relationships given by 
ds = a 2 cos^, sin ^d^dy and 0 = 7t - 2L >1 . On the other hand, according to the definition 
of scattering phase function denoted by P, one may define the following relationship: 


1 = 


4 nr' 


PI 0 > 


(15) 


where G s is the scattering cross section. Therefore, after the two polarization 
configurations are accounted for, the contribution of the external reflection to phase 
function may be obtained from comparison of Eqs.(14) and (15): 


2 a s 

HI*/, I 2 +l*i I 2 ) 

2 a-a ’ 


(16) 


where Q e and Q a are the extinction efficiency and absorption efficiency, respectively. 
The value of is 2 since the particle is very large whereas Q a is determined by the ratio 
of refracted energy to the incident energy associated with the projected area. Note that the 
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divergence factor introduced by van de Hulst 25 is absent in Eq.(16) because the rays that 
transmit through the sphere are ignored here. 

An improved diffraction formulation has been given by Yang and Liou 26 that, 
unlike the conventional method that is limited to scattering angle smaller than 90°, can be 
applied the entire scattering angular region. When the improved method is applied to the 
diffraction by a sphere, the scattering amplitude matrix 25 is given by 


% M 

^1 j 


{kay 2Ji(£asin0) 
4 £asin0 


COS0(1 + COS0) 
0 


0 > 1 


1 + COS0 


(17) 


where J, is the Bessel function of the first kind. The phase function that includes the 
contribution from diffraction and external reflection can be obtained from Eqs.(16) and 
(17) as follows: 


P(0) = 


n2 


(ka) 2 ’2./j(fcasin0) 
8(Qe~Qa)l ka sin 0 J 


(1 + cos 2 6)(1 + cos0) 2 + 1 (l ^ l - ^ 1 } • (18) 


2 Qe-Qa 

In addition to the phase function, the polarization configuration should be used to identify 
the correct combination of reflection coefficients. When the incident radiation is 
unpolarized, the polarization feature of scattered radiation can be specified by the degree 
of linear polarization (DLP). When scattered radiation is composed of diffraction and 
external reflection, the DLP associated with the phase function in Eq.(18) may be given 
as follows: 


DLP(0) = 


/±(Q) ~ -(v(0) + f±(9) ~ 1 // (®) 
/'(0) + #(0) + 71(0) + 7^(0) 


(to) 2 [2J 1 (fasin0)/tosin0] 2 (l + cos0) 2 (l-cos 2 0) + 4(l/? x l 2 -[ R„ | 2 ) 
(fo2) 2 [2/,(fcjsin0)/A:asin0] 2 (l + cos0) 2 (l + cos 2 0) + 4(|7? ± | 2 +| R// 1 2 )’ 
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if and I[ (if and If) are the scattered intensity associated with diffraction and external 
reflection, respectively, for the case that the electric field is polarized along a direction 
that is perpendicular (parallel) to the scattering plane. 

Figure 3 shows the phase function and DLP computed from Mie theory and the 
geometric optics method at a wavelength of 12 pm. The Mie computational code 
developed by Wiscombe 27 is used in this study. Since the size parameter is extremely 
large (%= 1047.2), the scattered energy at a large scattering angle (approximately 20° in 
this case) is primarily from the reflection of incident rays. Through comparison of the 
geometric optics solution to the Mie results, we can identify the correct reflection 
coefficient pair. Evidently, the reflection coefficient pair (R TM j,,R T e,±) produces a phase 
function and DLP consistent with the Mie solution at large scattering angels where the 
external reflection dominates. Since the Mie solution for DLP is positive for the entire 
scattering domain, the reflection of incident wave with vertical polarization configuration 
is much larger than the reflection for the case when the polarization is parallel to 
scattering plane. For this reason, if we incorrectly apply the reflection coefficient R m L 
for the incident wave polarized perpendicularly with respect to scattering plane, 
substantial errors can be produced for both phase function and DLP, as is evident from 
Fig.3. 

Figure 4 is similar to Fig.3, except that the wavelength is 1 1 pm and the real part 
of the refractive index is close to unity (see Fig. 1). It can been seen that the phase 
function values at large scattering angles at 11 pm are much smaller than the 
corresponding results at 12 pm since external reflection is weaker for the former 
wavelength. However, the correctness of reflection coefficients used in the computation 
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is also critical to the phase function and DLP even in the optically tenuous case. As a 
matter of fact, the errors associated when the wrong wave mode is assumed for the 
refracted wave are even more pronounced for the results shown in Fig.4 than for those 
shown in Fig.3. From these results, it is clear that the correct wave mode must be used 
for the refracted wave in the derivation of the reflection coefficients regardless of the 
magnitude of the real part of refractive index whenever the particle is strongly absorptive. 

The comparison of the ray-tracing method to Mie theory for the phase function of 
spheres has been presented by Liou and Hansen 28 for a polydispersive case having an 
assumed particle size distribution. For a monodispersive case (a system composed of 
particles having a single size and shape), the comparison study has been presented by 
Macke et al. 29 for spheroids and finite circular cylinders at a near infrared wavelength, 
which shows that the applicable size parameter for the geometric optics method is 
approximately 60. In addition, these authors pointed out that the size parameter required 
for a convergence of the exact solution and the geometric optics solution is larger for 
spheres than for nonspherical particles. 

The results of this study suggest that the geometric optics method can be applied 
to much smaller size parameters when strong absorption is involved. Figure 5 shows the 
phase function computed from Eq.(18) in comparison with the Mie solutions. Evidently, 
when the size parameter is larger than approximately 30, the geometric optics solution 
can be used to approximate the exact solution. For size parameters larger than 50, the 
geometric optics solution essentially converges to the exact solution in the side and 
backscattering directions. 
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From the preceding discussions associated Figs. 3, 4 and 5, it is clear that TM and 
TE modes should be applied to the components of electric field that are parallel and 
perpendicular to the incident plane, respectively. By using this combination of wave 
modes, it is straightforward to derive the refracted wave on the basis of electromagnetic 
boundary conditions. Through application of the ( TM„,TE ’ ± ) wave mode to the refracted 

wave, we can obtain the refracted electric field as follows: 


where 


E ,,o = E t,o,A + E ,,o$ + E ',o,£' 


T 0 

a ( E . ^ 

0 T 

u h e 

rrt Q V i0 $ ) 

1 V J 


T 2(AA r + iN i )cosC li 

a e cos + ( N r cos + iN n ) ’ 

T 2cosC,. 

p cos ^ ( . + (N r cos +iN n )’ 


i2N n sin^ cos £,• 
Ecos^, + (N r cos L, t + iN n ) 


(20a) 

(20b) 

(20c) 

(20d) 

(20e) 


Note that in Eqs.(20a)-(20e) the dependence of these parameters on the wave mode is not 
specified explicitly. It should be pointed out that the transmission matrix in the 
conventional ray-tracing scheme is a 2x2 matrix. When absorption is involved, the 
transmission matrix becomes a 3 x 2 matrix. When absorption of the medium is absent, 
T y is reduced to zero because N n is zero and Eqs.(20a)-(20e) are reduced to the 


conventional Fresnel refraction coefficients. 
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3. Asymptotic Solution for Scattering and Absorption Properties of Large 
Polyhedral Particles with Strong Absorption 

A. Physical Basis for Asymptotic Solution 

In the previous section, it is shown that the conventional geometric optics 

approach in terms of the superposition of diffraction and external reflection can well 

approximate the Mie solution for the spherical particle case when the size parameter is 

large and the particle is absorptive. Since the variation in the surface normal direction is 

continuous for a sphere, a continuous distribution of the externally reflected energy is 

obtained. However, this is not true if the particle geometry is a polyhedron such as a 

hexagonal plate or column, which is the basic geometric structure of ice crystals in cirrus 

clouds. In the numerical solution, there is a physically incorrect discontinuous 

distribution of scattered energy for any specific orientation of the particle. 

For example, consider the scattering of radiation by a hexagonal ice column when 

the incident radiation is perpendicular to the c-axis of the particle and parallel to a line 

that connects two symmetric apexes of particle cross section. If the rays refracted into the 

particle are ignored due to absorption, the phase function can be written as follows: 

sin[(£Z, / 2) sin 0 5 cos <p s ] sin[( V3Aa / 2) sin 0 5 sin (p J 
(&Z,/2)sin0 5 cos(p 5 (V3foa/2)sin0 s sincp 5 

(1 +cos 2 e j )(l +COS0J 2 + '^ mL [6(<p, - ft / 2) + 5(<p, + K / 2)] 

SO* -2*/3)](| V//.I I 2 + 1 S re.±,C,-itrt I 2 ). (21=) 

a, = 1 f - 1 Rr e . u ,„, t I 2 ). (21b) 

where L and a are the length and semi-width of hexagonal columns and 8(x) is the Dirac 
delta function. To avoid any possible confusion, in Eqs.(21a) and (21b) we use 0 5 and (p 5 


ne„<P,) = 


3 k^La 
8 na c 
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to indicate the scattering angle and the azimuthal angle of the scattering plane. Here the 
plane of (p 5 = 0° is the plane that contains incident direction and the c-axis of the 
particle. On the right side of Eq.(21a), the first term represents the contribution from 
diffraction and the second term represents the contribution of external reflection. 
Obviously, the second term is nonzero only at two angles as is expressed by the delta 
functions. This singularity is inherent in the conventional geometric optics method. The 
delta-reflection involved in Eq.(21) originates with the same physical mechanism of the 
forward delta-transmission that has been well explained by Takano and Liou 12 and 
Mishchenko and Macke 30 . It should be pointed out that the singularity in Eq.(21b) could 
be overcome numerically in the ray-tracing calculation if the particles are randomly 
oriented in space. The average intensity is given by the scattered energy divided by the 
corresponding solid angle elements. The discontinuity may be avoided given proper 
resolution for the solid angle elements with sufficient orientations of the particles, i.e., 
with sufficiently random orientation. 

To circumvent this particular disadvantage of the conventional method for a 
specific orientation, we apply the geometric optics method to solve for the internal field 
inside the particle following Yang and Liou 31 . The scattered far-field, extinction and 
absorption cross sections are given by 

E s (f) = - 1)JJJ ) - r[r.E(F )]} exp (-/Af.r)rfV , (22a) 


a ext = Im l 


|-(£^)Jp(F)£,(F) 


E: 


wV. 


® abs 


i I v 


(22b) 

(22c) 
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The integrals involved in preceding equations may be calculated using the ray-by-ray 
integration (RBRI) algorithm 31 . From the RBRI calculation and Eqs. (22a)-(22c), the 
amplitude scattering matrix and extinction and absorption cross sections are given by: 



direction; a s = e s xp f ; a o = e, x j3 f ; and /(r 1 ) is the ray penetration depth with the 

incident point at r ' . The other variables involved are defined in the previous discussion. 
If the absorbing particle is large, it is expected that exp[-£A(./(r' )] — > 0 . In this case, 

analytical expressions may be derived for the integration involved in Eqs. (23a) and (23c) 
if the surface of the particle is locally flat such as in the case of polyhedral geometry. 

For a given flat face of the particle surface, Fresnel reflection does not vary with 
the location at the face because the incident angles are the same for all the points on this 
face. Thus, the following matrix is independent of the position of a point on a flat surface: 

■ - , - Jz o 

S 2 S 3 ) cos£, a, -a, a, P 0 T 

~s< sj. Nr+iN, a, V s e s ) T Q p 


A \ 

“'f* . (24) 

P«. PP.J 
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Therefore, the integral in Eq.(23a) should actually be evaluated only with respect to the 
phase variation over the local flat faces of the particle. For evaluation purposes, the 
surface of a polyhedral particle can be divided into a number of area elements shaped as 
parallelograms and triangles. For a parallelogram or a triangular face, the position vector 
of a point on this face is expressed by 


r = 


r 0 + ^ + r|r 2 for parallelogram 
r o + E,(l - r|)^ + r)£,r 2 for triangle’ 


(25) 


where r 0 is the position vector of the apex of the geometric shape, and r ] and r 2 constitutes 
the neighboring two sides of the shape. The vector ^e[0,l] and r|e[0,l]. From 
Eqs.(24) and (25), the integration of phase variation over a parallelogram area element 
yields 


D< = JJ expire, ~e s )-rj}d 2 rj 


face j 


= l 0, x 0,2 1 f f ex P[^(0 “ «,) • (o.o + no, + 

Jo Jo 

=1 0, x 0.2 1 ex P Wfi, - e s ) • (o.o + 0, 1 2 + 0,2 1 2 )1 

sin[fc(e. - e 5 ) • Q., /2] sin[ft(e,. - ej • Q , /2] 

*(O-4)-0, /2 Ke,-e s )-rjf2 

Similarly, if the phase variation is integrated over a triangular area, it follows that 
D i = JJ exp{/£(e, -e s )- f^f } 


(26a) 


j 


exp[ifc(e,.-e>Q. 0 ] 

0,1 " 0,2 4) -(0,2-0,) 


=1 r. , x r, , I _ : ' <1 exp[/A:(e / -e s )r J2 / 2] 
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(26b) 


-exp[ik(e : -e s )-r jA /2] 


sin[A:(g. /2] | 

k{ei-e s )? h J2 


The scattering matrix that includes the contribution from all the faces can be written as 


S 2 (e s ) S,(e s f 

Me,) J 


^(e - 1) 

4 % 



r S 2 (e s ) 

U(4) $&)J/ 


(27) 


where h{e t ■ hj ) is a step function that indicates whether the face is illuminated by 


incident radiation: 


h{e r hj) = 


fl, for e i h j > 0 
{ 0, for e i ■ hj < 0 


(28) 


Once the amplitude scattering matrix is given, the calculation of phase function is 
straightforward on the basis of the definition presented by van de Hulst 25 . It should be 
pointed out that a special form of Eq.(27) have been given for hexagonal particle without 
an incorporation of the inhomogeneity effect of the refracted wave on the transmission 
matrix 31 . The result presented in this study is more general and can be applied to any 
polyhedron with strong absorption. 

Since the multiple scattering calculation is based on the solution of the radiative 
transfer equation, the most important parameter of concern is the single-scattering albedo 
rather than extinction or absorption cross section. In particular, in remote sensing 
applications, a pre-described optical thickness is usually used in the calculation of look- 
up tables for the microphysical or optical properties of cloud particles. For this reason, in 
the present calculation we focus on the phase function and single-scattering albedo. The 
later parameter is defined as 


® ext ® a bs 

®ext 


(29) 
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where the extinction and absorption cross sections, <5^ and G abs can be calculated from 
Eqs. (23b) and (23c). 

B. Numerical Results for Hexagonal Ice Plates and Columns 

As an application for the asymptotic solution for light scattering by a strongly 
absorptive particle, we present numerical results for ice crystals composed of hexagonal 
plates and solid columns. A hexagonal shape is selected here because this shape has been 
studied extensively for inference of cloud optical thickness and particle size in remote 
sensing applications (e.g., Minnis et al. 32 , Baum et al. 33 , Han et al. 34 ). 

Fig.6 shows the comparison of the phase function computed by the FDTD method 
and the asymptotic theory (Eq. 27) for hexagonal ice crystals of two sizes: 2a/L= 25pm 
/25pm and 40pm/40pm, where L and a are the length and semi-width of the hexagonal 
crystal. A random orientation condition is assumed for ice crystals in both the asymptotic 
and the FDTD calculations. It can be seen from Fig. 6 that the asymptotic solution more 
closely approximates the FDTD results for the larger particle size of 40/40 pm. While 
one might infer that the asymptotic theory can be applied to hexagonal ice crystals with 
sizes larger than approximately 40pm in the infrared window region, the convergence 
between FDTD and asymptotic theory will, in fact, depend on magnitude of absorption. 
For large particles at IR wavelengths, the FDTD calculation is computationally 
expensive. Thus, the asymptotic theory provides an efficient way to estimate the phase 
function of large particle with strong absorption. This advantage of the asymptotic theory 
in the calculation of phase function may be useful in practice. For example, the 
parameterization of cloud longwave radiative properties requires the phase function 
information though the use of either an asymmetry factor or a scaling factor 35 , which can 
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be provided by the asymptotic theory efficiently and accurately when the particle size is 
larger than approximately 40 pm. 

Figure 7 shows a comparison of the single-scattering albedo computed from the 
FDTD and the asymptotic theory. In this diagram, the results computed for the equivalent 
spheres are also shown. Following Mitchell et al. 36 , Fu et al. 37 , and Grenfell and Warren 39 , 
we define the radius of the equivalent sphere for a hexagonal particle as follows: 

3 V 3 V3/2 aL 

r e = = — 7= , (30) 

e 4 A 2V3/2 a + L 

where V is the volume of nonspherical particles, and A is the projected area. Since strong 
absorption is assumed in the computation of the single-scattering albedo in the 
asymptotic theory, the result does not change with the variation of particle size. As a 
result, the asymptotic theory can lead to substantial errors in single-scattering albedo 
calculation for small and moderate particle sizes. However, from the trend of the FDTD 
solution, it is expected that a convergence between the FDTD and asymptotic results can 
be reached when the particle size is on the order of 200. Approaches developed by Fu et 
al. 37 , Michell 39 , and Baran et al 40 can be used to overcome the shortcomings of the 
asymptotic theory in calculating single-scattering albedo. It should be pointed out that 
these approaches cannot provide the information for the phase function. For this reason, a 
combination of these approaches (for the computation of single-scattering albedo) with 
the asymptotic theory (for the computation of phase function) may be more useful in 
practice. 

Figures 8 and 9 show the phase functions of randomly oriented ice columns and 
plates in comparison with the phase function of equivalent spheres. At large scattering 
angles, the phase function values are essentially the same for the spherical and 
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nonspherical particles. The external reflection becomes insensitive to particle geometry 
once the random orientation condition is assumed. However, for moderate particle sizes 
(r e =26J pm) , substantial differences may be noted for the phase function at scattering 
angles in the 1 0°-60° region for the two particle shapes. 

Figure 10 shows the variation of the single-scattering albedo of hexagonal ice 
plates and columns that are randomly oriented in the infrared region between 8 pm and 
16 pm. Also shown are the real and imaginary parts of refractive index data compiled by 
Warren 41 . Evidently, the variational pattern of the single-scattering albedo follows that of 
the real part of the refractive index. As we have mentioned, the real part of refractive 
index has a value close to unity in the vicinity of 1 1 pm, and the corresponding imaginary 
part is substantially large. Near 1 1 pm, the extinction of incident radiation is dominated 
by absorption, as is known as the Christiansen effect 23 . As a result, a minimum is noted in 
the single- scattering albedo. From Fig. 10, it can also be seen that the single-scattering 
albedo for hexagonal particles is larger than for spheres. This feature has been noted in 
the comparison of the FDTD solution with the equivalent spherical results that are shown 
in Fig. 7. 

It is a common practice to assume ice crystals are randomly oriented in space. 
However, in the atmosphere specific orientations are preferred for large plates and 
columns. For plates in an environment of a low Reynolds number, such as in a typical 
midlatitude cirrus cloud, the a-axis (the hexagonal asymmetric axis) of the particle tends 
to face vertically whereas the oaxis tends to be oriented horizontally. In other words, the 
crystals tend to align themselves so as to reach a stable state in its falling process. At 
visible wavelengths, the effect of preferred particle orientation on the phase function has 
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been discussed by Takano and Liou 42 and Rockwitz 43 . As pointed out in Section 3. A, the 
application of a conventional ray tracing technique will produce a discontinuity of the 
scattered energy and consequently cannot produce a continuous phase function. This 
shortcoming is avoided with the present analytical asymptotic theory. When preferred 
orientations are assumed for ice crystals, the phase function depends not only on the 
scattering angle but also on the scattering azimuthal angle, i.e., on which scattering plane 
the scattered field is observed. Figure 1 1 shows the geometry for the scattering 
configuration related to a plate oriented with its c-axis aligned vertically. The principal 
plane ((p s = 0°) is defined as the plane that contains the incident direction and zenith 
directions. The plane that is normal to the incident direction is denoted by the dotted 
ellipse whereas an arbitrary scattering plane is denoted by the dashed ellipse in the 
diagram. On a given scattering plane, the positive scattering angle is measured, with a 
view along the incident direction, clockwise from the principal plane. For example, on 
the principal plane the scattering angle is positive if it is measured from forward 
scattering direction to the vertical zenith whereas the scattering angle is negative if it is 
measured from forward scattering direction to nadir view. 

Figure 12 shows the phase function for hexagonal plates having a preferred 
orientation for four scattering azimuthal angles. Ice crystals are rotated randomly about 
their c-axis. The strong forward peak caused by diffraction can been seen for the four 
cases with different scattering azimuthal angles. For the case of <p 5 =0°, very strong 
scattering peaks may be noted at scattering angle of 120° and 60° for incident zenith 
angles of 30° and 60°, respectively. The scattering maxima at 120° and 60° correspond to 
the specular reflection from the top faces of the plates. For cp 5 = 0°, 30° and 60°, 
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significant asymmetry in the phase function is observed in the scattering angle regions of 
(0°, 180°) and (-0°, -180°). However, symmetry is observed for cp 5 = 90°, as is expected. 
For plates with a preferred orientation, there are a number of ripple structures in the phase 
function that are caused by phase interference. 

Figure 13 shows the phase function for ice columns with preferred orientations. 
The columns, randomly rotating about the zenith direction, are horizontally oriented but 
rotate randomly about their c-axis. The ripple structure noted in Fig. 12 for plates due to 
the phase interference is largely smoothed out in the case of oriented columns because the 
orientation of the particles in the column case is more random. Similarly to the case in 
plates, the specular reflection feature is apparent in the case where (p 5 = 0°. 

Figure 14 shows the single-scattering albedo corresponding to the phase functions 
in Figs. 12 and 13. Compared with the results for randomly oriented particles, the single 
scattering for both oriented plates and columns are smaller if the incident zenith angle is 
less than approximately than 60°. For plates, there is a pronounced peak in the single- 
scattering albedo at the incident zenith of 70°. This is thought to be due to the external 
reflection for the top face of the plate increasing with an increase of incident zenith angle. 
However, when this angle approaches 90°, the projected area of the top face is essentially 
zero and the external reflection from the side faces, along with diffraction, dominates the 
scattered field. 

4. Summary and Conclusions 

An asymptotic solution in an analytical format is presented for the scattering 
properties of a general polyhedral particle with strong absorption. For scattering 
calculations involving nonspherical particles such as aerosols and ice crystals, we show 
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that the transverse wave condition is not applicable to the refracted electromagnetic wave 
when absorption is involved. In the geometric optics solution for wave propagation in an 
absorbing dielectric medium, either the TM wave condition (i.e., where the magnetic 
field of the refracted wave is transverse with respect to the wave direction) or the TE 
wave condition (i.e., where the electric field is transverse with respect to the propagating 
direction of the wave) may be assumed for the refracted wave to satisfy the 
electromagnetic boundary condition. The wave mode assumed for the refracted wave 
affects both the reflection and refraction coefficients. As a result, a nonunique solution 
for these coefficients is derived from the electromagnetic boundary condition. Through 
comparison of the asymptotic geometric optics solution to results determined from Mie 
theory for absorbing spheres, it is shown that TM and TE wave modes should be applied 
to two polarized components that are parallel and perpendicular to the incident plane, 
respectively. In this study we have identified the appropriate solution for the Fresnel 
reflection coefficients in the ray tracing calculation. We present the 3x2refraction or 
transmission matrix that completely accounts for the inhomogeneity of the refracted wave 
in an absorbing medium. Based on the Fresnel coefficients for an absorbing medium, we 
derive the asymptotic solution in an analytical format for the scattering properties of a 
general polyhedral particle. Numerical results are presented for hexagonal plates and 
columns with both preferred and random orientation. The asymptotic theory produces 
reasonable accuracy in the phase function calculations in the infrared window region 
(wavelengths near 10 pm) if the particle size is on the order of 40 pm in diameter or 
larger. 
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However, since strong absorption is assumed in the computation of the single- 
scattering albedo in the asymptotic theory, the single scattering albedo does not change 
with variation of the particle size. As a result, the asymptotic theory can lead to 
substantial errors in the computation of single-scattering albedo for small and moderate 
particle sizes. However, from comparison of results with the FDTD solution, it is 
expected that a convergence between the FDTD results and the asymptotic theory results 
can be reached when the particle size approaches 200 pm. 

For two infrared wavelengths at 1 1 and 12 pm, we show that the phase function at 
side- scattering and backscattering angles is insensitive to particle shape if the random 
orientation is assumed. For preferred orientations, however, we show that the phase 
functions for plates and columns are significantly different. Additionally, when a 
preferred two-dimensional orientation condition is assumed, the phase function is 
observed to have a strong dependence on scattering azimuthal angle. Moreover, 
numerical results show that the single-scattering albedo has a very strong dependence on 
the inclination angle of incident radiation with respect to the rotating axis for the 
preferred particle orientations. 


This research has been supported by a grant of NASA’s MODIS project and 
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Figure Captions 

Fig.l. Comparison of the inherent and apparent refractive indices for ice at 1 1 and 12 pm. 

Fig.2 Scattering geometry for a large sphere with strong absorption. 

Fig.3 Phase function and the degree of linear polarization computed from geometric 
optics in comparison with Mie solution at X=\2 pm. 

Fig.4 As the same as Fig.3 except that the wavelength is 1 1 pm in Fig.4 

Fig.5. Comparison of the phase functions computed from Mie theory and the geometric 
optics method for three moderate sizes at wavelengths 1 1 and 12 pm. 

Fig. 6 Comparison of phase function computed from the FDTD technique and the present 
asymptotic theory for hexagonal ice crystals. 

Fig. 7 Comparison of the single-scattering albedo computed from the FDTD technique 
and the present asymptotic theory for hexagonal ice crystals. Also shown are the 
results computed for equivalent spheres. 

Fig. 8 Comparison of phase function of hexagonal columns with that of equivalent 
spheres. 

Fig.9 As for Fig.8, except for hexagonal plates. 

Fig. 10 The comparison of the single-scattering albedo values computed for hexagonal 
ice crystals and equivalent spheres. 

Fig. 1 1 Incident and scattering geometry for an ice crystal with a preferred orientation. 

Fig. 12 Phase function computed for hexagonal plates with preferred orientations. 

Fig. 13 As for Fig. 12 except for columns. 

Fig. 14. The single-scattering albedo for ice crystals that have preferred orientations in 
comparison with the results for the randomly oriented case. 


36 


Incident Angle (deg.) Incident Angle (deg.) 


Real Part of Refractive Index 



Inherent Refractive Index : (m r +imj), Apparent Refractive Index: (N r +iNj) 











A,=12 |im, a=2000 jim (ka =1047.2), m r +inrij=1.2798+i0.41333 



Fig.3 






Phase Function DLP Phase Function 


>*=1 1 jLim, a=2000 (ka =1 142.4), m r +irrij=1.0925+i0.248 


Mie Theory 

~ ]q 6 -| ( r te,//»Rte,i) 

c Mode 


c 

10 7 

o 

V- ■ 
o 

10 6 

a 

3 

10 5 

UL 

10 4 

a> 

c/> 

10 3 

co 

-C 

10 2 

CL 

10 1 


0 1 2 3 4 5 

Scattering Angle ( 0 ) 


8 io _1 


Mie Theory 


( r tm,/ 

Mode 


//> r tm,±) 


I io j' 

io 2 - 

.CL 1Q 1 - 


0 1 2 3 4 5 

Scattering Angle ( 0 ) 













120 180 0 60 120 180 

Scattering Angle ( deg.) 


T1 

(Q 

Ol 


Phase Function 


q 

rb 


O 

o 


o 

ro 


o 

co 



o 

rb 


o 

o 




O 

CO 




O 

i 

ro 


O 

o 




o 

CO 


s. 



o 

« 

ro 


o 

o 








Phase Function Phase Function 











Func 












>2 






Wavelength ( nm ) Wavelength ( nm ) 




Single-Scattering Albedo 



<y> 







W'BU 



/ "S ,.N V 

^ / V 


g o 13 

-T Q ) 

CD 3 3 
2- Q. CD 

O S’ 3 
3 CD O 

—5 — \ 

~ 3 

D) 


. CD 



0 ) “0 
O 3. 
CD 3 
!=t O 

CD - o ' 


- nC 0 > 

Isa 

3 !=fc =*■ 
CD CD £D 

5'^ 

CQ 


Incident Rays 







120 -60 0 60 120 180 -180 -120 -60 0 60 12 

Scattering Angle (deg.) Scattering Angle (deg.) 



Incident Radiation 














